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Two computer programs, known as Matrix Through-Flow and Matrix 
Blade-To-Blade, for analyzing the meridional and blade-to-blade flow 
patterns are described. The numerical solutions are obtained by finite 
difference approximations to the governing Poisson-type differential 
equations for the stream function. Solutions for several turbomaehines, 
giving flow patterns and velocity distributions, are included. 


The flow through a modern gas turbine or compressor is an extremely 
complicated three-dimensional phenomenon. The flow has strong gradients 
in the three physical dimensions — axial, radial, and circumferential — as 
well as time and viscosity effects. The observation that the flow problem 
was not easily amenable to numerical solution led early investigators to 
search for a design system having ease of application. The computational 
difficulties were resolved by making approximations which permitted the 
use of two-dimensional techniques. These approximations were based on 
two flow models, 

(1) Blade element flow 

(2) Axially symmetric flow. 

The blade element approach assumes that the flow in the blade-to-blade 
or circumferential plane can be described by considering the flow around 
blade profiles formed by the intersection of a cylindrical flow surface and 
the blading. 

Axial symmetry assumes that an average value can be utilized to 
represent the state of the fluid in the blade-to-blade plane. 

On the basis of these two flow models, several investigators developed 
analysis and design methods for the axial-flow compressor and turbine. 
In the case of the compressor, one of the earliest design methods appeared 
in Howell’s classic papers in 1945 (refs. 1 and 2) . Using the blade element 
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flow model, Howell correlated experimental linear cascade data to es- 
tablish a limit that has to be placed on the allowable deflection in any one 
blade row and determined empirical rules for the deviation and flow loss. 
In estimating the overall performance of the compressor, the flow is 
analyzed along a “mean” or “reference” diameter and the gas state is 
estimated at planes between adjacent blade rows, making use of the 
axially symmetric flow model. Similar methods were developed for the 
axial-flow turbine and, of these, the method of Ainley and Mathieson 
(ref. 3) is one of the best known. These relatively simple, albeit one- 
dimensional methods for analyzing the overall properties of the flow field, 
developed when the digital computer was in its infancy and the develop- 
ment of methods suitable for hand desk machines was one of the prime 
goals, are still, in principle, used widely throughout the aircraft industry 
and are likely to remain in use for some time. 

More recently, with the advent of the large, high-speed digital com- 
puter, techniques (refs. 4 and 5) have been developed for analyzing the 
subsonic fluid motion in the meridional or hub-to-tip plane of axial-flow 
machines at stations other than the mean diameter (which was used in 
the early days) both inside the blade rows and in the duct regions. Similar 
methods have been developed for centrifugal and mixed-flow impellers by 
Hodskinson (ref. 6) and Wood, et al. (ref. 7). In parallel, several in- 
vestigators (refs. 8, 9, and 10) have been working on the problem of 
generating a computer solution for the subsonic blade-to-blade flow with 
allowances for radial acceleration imposed by the curvature of the stream- 
lines in the meridional plane and for the effects of Coriolis forces. 

The purpose of this paper is to present an outline of tw'o advanced 
computer solutions that have been developed at the National Gas Turbine 
Establishment (NGTE) for the meridional and blade-to-blade flow 
patterns. Solutions for several turbomachines, giving flow patterns and 
velocity distributions, are included. 

MATHEMATICAL ANALYSIS 

The mathematical analysis is based on the earlier work of Wu (ref. 11) 
who developed a general theory for the three-dimensional, inviscid, 
steady flow through an arbitrary turbomachine. The equations of motion 
are satisfied on two intersecting families of stream surfaces known as the 
first kind, SI (blade-to-blade), and the second kind, S2 (meridional), the 
complete flow solution being obtained by an iterative process between the 
flows in the two stream surfaces. 

Stream Function Equation for SI Surface 

In the real blade-to-blade flow, the SI stream surface would be twisted. 
To permit computations of the potential flow in the blade-to-blade plane 
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of a stationary or rotating blade row, Smith (refs. 10 and 12) assumed 
the stream surface was a surface of revolution. 

The shape of the SI surface is obtained by rotating a streamline in the 
meridional plane (fig. 1) about the axis of rotation. 1 In order to analyze 
the flow through any type of turbomachine, it is convenient to rotate the 
r,z axes through an angle 6. Using x and as the two independent, vari- 
ables, the continuity equation and the equations of motion can be mani- 
pulated to arrive at a Poisson-type differential equation for the stream 
function. 
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where the stream function \p is defined by 
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and the velocity components W x and W v are related by 

W v = — W x tan X (3) 

Equation (3) is the geometrical condition that the flow follows the stream 
surface. The derivatives in equations (1) and (2) are those which Wu 



Figure 1. — Meridional plane. 



1 For two-dimensional cascade flow the stream surface is a cylinder. 
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refers to as special derivatives taken on the stream surface, and the 
integrating factor b is proportional to the local thickness of a thin stream 
sheet whose mean stream surface is the 51 surface considered. 


Stream Function Equation for S2 Surface 


To analyze the flow in the meridional plane of an arbitrary turbo- 
machine, Marsh (ref. 4) developed a matrix through-flow method. Inside 
the blade rows, the flow is analyzed in an 5 2 stream surface and for the 
duct regions between adjacent blade rows, the flow is assumed to be 
axially symmetric. 

As in the case of the 51 surface, the r,z axes (fig. 1) are rotated through 
an angle 6 and x,y are the two independent variables. In a manner similar 
to the 51 solution, an equation for the stream function can be derived. 
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where the stream function satisfies 
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The integrating factor B in equation (5) is proportional to the local angu- 
lar thickness of the 52 stream surface and in the through-flow analysis it 
is assumed to be proportional to the width of the blade passage. In formu- 
lating the stream function equation — equation (4) — the viscosity terms 
were omitted in the equations of motion but the entropy terms were 
included, and Marsh introduced the effects of irreversibility into the flow 
calculation by defining a local polytropic efficiency for expansion and 
compression. 

For the flow to follow the stream surface, within the blade rows, the 
three components of velocity are related by 


Wt= — IF* tan X — IF* tan p (6) 

In the duct regions there is no change of angular momentum along a 
streamline and the circumferential velocity satisfies the relationship 


rV* = constant 


(7) 
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NUMERICAL SOLUTION 

The equations for the stream function — equations (1) and (4) — are 
nonlinear, but they can be solved iteratively using finite difference 
techniques. 

Finite Difference Approximations 

In conventional finite difference analysis the domain is covered with a 
square or rectangular grid and a five-point star is used since this leads to 
a simple approximation for the Laplacian operator. However, for an- 
alyzing the flow in the SI and S2 stream surfaces such a simple grid is not 
accurate enough, owing to the irregular boundaries of the flow domain 
giving rise to boundary finite difference stars with short limbs and con- 
sequently a large truncation error. A good example of this, in fluid 
mechanics, is the recent blade-to-blade method developed by Katsanis 
(refs. 8 and 9) in which the flow domain is covered with a square grid. 
It is clear that the truncation error is significant since the boundary 
condition of zero velocity normal to the blade surfaces is not satisfied. 

In the NGTE methods, use is made of the powerful software of present- 
day digital computers by adopting an asymmetric finite difference grid. 
The grid (fig. 2) consists of straight lines normal to the x direction, each 
line having the same number of equally spaced grid points. In the case of 
the SI surface, the blade suction and pressure surfaces form curved grid 
lines, and for the S2 surface, the inner and outer annulus walls form 
curved grid lines so that there are no additional difficulties for grid points 
close to the boundaries. The spacing of the straight lines need not be 
uniform and where necessary can be varied locally (in the blade leading 
and trailing edges, for instance) in order to obtain a detailed picture of 
the flow. 



BLADE -TO -BLADE 
COORDINATES 


MERIDIONAL 

COORDINATES 


Figure 2.— Asymmetric grid. 
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To formulate a finite difference approximation for the Laplacian 
operator having an error of fc 2 , where fc is the local grid spacing, equations 
(1) and (4) were modified by adding the term E (d\p/dx) to both sides. 
Thus, for the Si surface the stream function equation becomes 


dx 2 r 2 d<jr 


d\p 

dx 


= F(x,4>) -\-E 


d\p 

dx 


or 

d\I/ 

vy+E^ = q{x,<l>) ( 8 ) 

dx 

where E is a function of the grid spacing in the x direction and is zero for 
uniform spacing. The operator Vty+E (dxp/dx) is approximated by a 
ten-point star for the interdependence of the function values at neigh- 
boring grid points. To maintain an overall accuracy of order fc 2 , the 
derivative d\p/dx is also approximated by the use of a ten-point star. 

Boundary Conditions 

Considering first the S2 surface, the boundary conditions are relatively 
simple. At inlet to the turbomachine, the flow conditions are known; 
therefore, the stream function distribution is defined for the first straight 
line of the grid. The inner and outer annulus walls form limiting stream- 
lines, so that for grid points on the walls the stream function is known. 
For the far downstream boundary, it is assumed that the shape of the 
exit duct is such that the stream function distribution is the same on the 
last two straight lines of the grid. 

The blade-to-blade problem — Si surface — poses quite complex bound- 
ary conditions. Far upstream of the blade row the gas state and flow angle 
are known and it is assumed that the flow is uniform. The gradient of 
stream function is defined, therefore, for the first straight line of the grid. 
Thus, from equation (2) 

Q 

= A <f> tan a u 

r « 

where A6 = 2ir/N and N is the number of blades in the row. For the blade 
region the suction and pressure surfaces form, by definition, limiting 
streamlines so that for grid points on the blade surfaces the stream 
function is known. Upstream and downstream of the blade the locations 
of the streamlines are not known until the problem is solved. For these 
regions, the boundary condition is that there is a circumferential perio- 
dicity of the flow. The final condition is that for the far downstream 
boundary. In a real blade-to-blade flow the circulation, and consequently 
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the outlet flow angle, is largely controlled by viscosity. In a potential flow 
model a criterion has to be adopted for fixing the circulation. In the 
method developed at NGTE, it is assumed that the flow is uniform and 
the flow angle is known far downstream of the blade row. These conditions 
fix the gradient of stream function on the last straight line of the grid 


r — „. „ /o'! 

vviiiuu, uuiu AO 



= — A<f> tan ad 
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Solution of Banded Equations and Convergence 

By making use of the finite difference approximations and the boundary 
conditions, the modified stream function equations — equation (8) for the 
(SI surface — can be written in matrix form: 

[M]-|>]=[g] (9) 

where [i p~\ and Q/] are column vectors formed by i and q at each grid 
point and [M] is a band matrix of the influence coefficients of the finite 
difference approximations. The method of solving equation (9) for the 
stream function is to solve for a given vector [g], to correct [g] using the 
new flow pattern, and then to repeat the cycle of calculation until the 
solution has converged to a specified tolerance. Since the matrix [M] is 
“banded,” only the band of nonzero elements is formed and stored in the 
computer and a very efficient direct method (ref. 13) is used to solve 
equation (9) for a given vector [g]. This method is better than the 
alternative indirect or relaxation method, as used by Katsanis, for the 
simple reason that it is very stable numerically. 

Numerical stability can be a major problem with any iterative method. 
In the matrix through-flow and blade-to-blade methods, the iterative 
process has been made stable by introducing a relaxation factor R ; thus, 

'Pp = 'Pp-i+R('l'-'f'p-i) (10) 

where 

calculated value for the pth iteration 
4 V value taken for the pth iteration 
4'p-i value taken for the (p — 1 ) th iteration. 

Additional stability was obtained in the through-flow method by limiting 
the percentage change in 4 / between successive iterations, a restriction 
which is automatically removed as the solution converges. For the blade- 
to-blade method, the stability was further improved for compressible flow 
by adopting a “marching” process of increasing the inlet Mach number 
gradually to the required value. 
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When the stream function is known, it is possible to calculate the 
products pW x , pW 4 , and pW y . To calculate the density and hence the 
velocity components a tabular method, as developed by Wu, is used in the 
through-flow method. For the blade-to-blade problem, an alternative 
method, suggested by Gelder (ref. 14), is used. In this method, the cal- 
culation of density is allowed to lag the stream function calculation by 
one iteration. This has the effect of improving stability and for com- 
pressible flow, the relaxation factor R — equation (10) — is a function of 
the maximum Mach number. 


BLADE-TO-BLADE FLOW PATTERNS 

Eight examples are given to illustrate the use of the blade-to-blade 
computer program. 

(1) Impulse turbine cascade 

(2) Seventy-degree camber blade 

(3) Axial turbine rotor tip section 

(4) Axial turbine rotor root section 

(5) Axial turbine stator blade 
(G) Turbine stator cascade 

(7) Three-dimensional flow past turbine stator blade 

(8) Radial cascade diffuser 

Impulse Turbine Cascade 

The first example is the incompressible flow past a 1 12-degree camber 
blade in cascade. The blade profile (fig. 3) is an impulse-type turbine 
blade having a pitch/chord ratio of 0.59 and 101 degrees flow deflection. 



Figure 3. — Blade profile-exact 
solutions. 


(b) 70 DEGREE CAMBER BLADE 
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The velocity ratios are plotted in figure 4. (Velocity ratio is defined as the 
ratio of local surface velocity to far downstream velocity.) Also shown is 
an exact solution obtained by Gostelow (ref. 15). The matrix solution is 
in very good agreement with the exact solution. 

Seventy-Degree Camber Biade 

This blade profile (fig. 3) has a pitch/chord ratio of 0.9 and 70 degrees 
of camber. The two-dimensional, incompressible velocity distribution for 
— 70 degrees of incidence is compared with an exact Gostelow solution in 
figure 5. In general, the matrix solution is in excellent agreement with the 
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PER CENT CHORD degree camber blade. 
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exact solution. It is noticeable that the main discrepancies are in the 
region of the blade trailing edge on the suction surface. For this region, 
the exact profile coordinates are a long way apart and it is probable that 
errors in interpolating the coordinates for the matrix solution have caused 
the discrepancies. There seems no reasonable doubt that complete agree- 
ment would have been obtained if the exact airfoil shape had been more 
fully defined. This example shows that there is no problem in analyzing 
high-incidence flows. The streamline pattern, calculated by the matrix 
method, is shown in figure 6. It may be seen that the leading edge stagna- 
tion point is well round on the suction surface. 

Axial Turbine Rotor Tip Section 

This example of two-dimensional, incompressible flow past a rotor tip 
section is given to illustrate the type of detailed flow pattern that can be 
calculated. The blade section is typical of a high pressure ratio turbine 
stage and is formed by a parabolic camber line and an analytical thickness 
distribution (ref. 16) . Initially, the profile was designed so that the blade 
inlet angle was equal to the gas inlet angle of 18 degrees, a condition often 
referred to as zero geometric incidence. Figure 7 shows the blade profile. 
The surface velocity distribution around the blade leading edge is plotted 
in figure 8. It may be seen that it has the undesirable characteristic of a 
high peak on the suction surface. Such effects have been found by Hall 
(ref. 17). This is due to the large induced incidence which can be seen 
from the streamline pattern in figure 9a. The high suction peak was 
reduced by effectively drooping the nose of the blade (fig. 7) by 10 degrees 
so that the profile was operating at — 10 degrees geometric incidence. The 
resulting velocity distribution is shown in figure 8 and, from the streamline 
pattern (fig. 9b), it may be seen that the induced incidence was consider- 
ably reduced. These results serve to show that computer methods can be 
very powerful in analyzing detailed aspects of the flow which would prob- 
ably be very difficult to find experimentally. 



Figure 6. — Streamline pattern for leading 
edge of 70-degree camber blade. 
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Axial Turbine Rotor Root Section 

This blade section is typical of a high pressure ratio turbine stage. The 
basis for the design is the same as that of the previous example. The gas 
inlet angle was 48.9 degrees and the blade geometry at inlet was chosen so 
that the geometric incidence was zero. At outlet, the blade passage was 
adjusted to satisfy the gas outlet angle of —63.9 degrees by the empirical 
rule of Ainley and Mathieson (ref. 3) . The blade surface velocity for two- 
dimensional, incompressible flow (fig. 10) shows that a detailed solution 
can be obtained in the region of the leading edge stagnation point. A 
particularly interesting feature of this blade section is that, according to 
the Ainley and Mathieson rule, the deviation 2 is 2.87 degrees negative. 
The velocity distribution for the trailing edge region is shown, enlarged, 
in figure 11 for an outlet flow angle of —64.15 degrees — a difference of 
only 0.25 degrees from the Ainley and Mathieson value. It is seen that on 
both the suction and pressure surfaces there is a rapid rise in velocity as 
the flow passes around the trailing edge. At the blade cutoff points, the 
velocities are equal, a criterion often used for fixing the outlet flow angle 
(ref. 18) . Also, if the two surface velocity distributions are extrapolated 
then the loading at the blade trailing edge is zero, thus satisfying Preston’s 
theorem (ref. 19) that equal and opposite vorticity should be shed from 


Figure 10 . — Velocity distribution 

turbine rotor root section. 



1 Deviation is the difference between the fluid and blade outlet angles. 
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Figure 11. — Trailing edge velocity 
distribution for turbine rotor root 
section. 


the two blade surfaces. The streamline pattern is shown in figure 12. It 
is seen that the flow leaves the trailing edge smoothly. This example 
shows that, applying existing velocity distribution criteria, the potential 
flow model gives an outlet flow angle in good agreement with well- 
established empirical rules, although it is perhaps surprising to find that 
the deviation is negative. 

Axial Turbine Stator Blade 

This example of two-dimensional, compressible flow is for the mean 
diameter section of a stator for a NASA turbine (ref. 20) operating at the 
design mass flow. The theoretical and experimental distributions of blade 
surface Mach number are compared in figure 13. In general, the computed 
Mach numbers agree well with experimental data. As mentioned earlier, 
Katsanis has developed a similar blade-to-blade method. In his recent 
paper (ref. 21) mention was made of an attempt to analyze the flow past 
this blade. He found that it was not possible to obtain an exact solution 3 
and he had to resort to an approximate solution. 

Turbine Stator Cascade 

This turbine cascade was fitted with blades having the same profile as 
the mean diameter section of the second-stage stator blades of the turbine 


3 The term “exact” has been used as meaning a numerical solution from the com- 
puter program. 




COMPUTER SOLUTIONS OF WU’S EQUATIONS FOR COMPRESSIBLE FLOW 57 


• EXPERIMENT FIGURE 

MATRIX 


13. — Mach number distribution for 
NASA axial-flow turbine. 



Figure 14. — Turbine cascade. 



Three-Dimensional Flow Past Turbine Stator Blade 

The flow through the two-stage turbine mentioned in the previous 
example has been analyzed using both the blade-to-blade and through- 
flow programs. The results of the through-flow analysis are presented in a 
later section of this paper. 

Two matrix blade-to-blade solutions for the flow past the second-stage 
stator blades were computed. The first solution was for two-dimensional 
flow (i.e., cylindrical stream surface of constant thickness) and the outlet 
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flow angle was calculated from the Ainley and Mathieson empirical rule. 
The second or quasi-three-dimensional solution is a refinement in that the 
stream surface thickness was varied. The variation of thickness was 
determined from a solution for a meridional flow pattern using the 
through-flow program and the outlet flow angle was determined by 
applying the condition of zero trailing edge loading. 

A comparison of observed blade surface Mach numbers with the 
theoretical calculations, for the mean diameter section, is shown in figure 
15. The most striking point here is that when some of the interactions 
between the meridional and blade-to-blade flow patterns are introduced 
the quasi-three-dimensional solution is in good agreement with experi- 
mental data. As mentioned earlier, this mean diameter section has been 
tested in cascade. The cascade Mach number distribution shown in figure 
14a corresponds to the turbine flow conditions given in figure 15. By 
comparing the cascade and turbine results, it may be seen that the three- 
dimensional flow effects are significant on the peak surface Mach number. 

Radial Cascade Diffuser 

To illustrate the types of turbomachines to which the matrix blade-to- 
blade method can be applied, the last example is a radial cascade diffuser. 
The initial calculations were made for incompressible flow with the cascade 
operating at zero geometric incidence and the outlet flow angle equal to 
the blade outlet angle (i.e., zero deviation). The theoretical velocity 
distribution is shown in figure 16. The peak near the trailing edge is due to 
the potential flow model picking up the rapid change in blade surface 
curvature in this region. In real flow, such peak velocities would be 
removed by the presence of boundary layers. By extrapolating the suction 

• EXPERIMENT 

TWO DIMENSIONAL ) MATRIX 

THREE DIMENSIONAL J SOLUTIONS 
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Figure 16 . — Velocity distribution for a 
radial cascade diffuser. 


surface velocity distribution from 90 percent of the way along the surface 
(i.e., just upstream of the rapid acceleration in the trailing edge region), 
it is seen that the condition of zero trailing edge loading is not satisfied. 
By increasing the outlet flow angle such that the deviation is 4.5 degrees, 
it is seen that the loading at the trailing edge satisfies Preston’s theorem. 

MERIDIONAL FLOW PATTERNS 

In this section, four examples are given of meridional flow patterns 
obtained from the matrix through-flow program. 

( 1 ) Two-stage axial-flow turbine 

(2) Single-stage axial-flow turbine 

(3) Low pressure ratio centrifugal compressor 

(4) High pressure ratio centrifugal compressor 

Two-Stage Axial Flow Turbine 

This turbine (ref. 22) is the one referred to in the previous section. In 
applying the matrix through-flow program, the effects of irreversibility 
were taken into account by assuming that the local polytropic efficiencies 
were constant throughout the flow field. From the comparison of the 
experimental and predicted profiles of axial velocity at the turbine exit 
shown in figure 17, it is seen that the through-flow theory gives a fair 
estimate of the axial velocities. Recent work by Gregory-Smith (ref. 23) 
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Figure 17 . — Axial velocity profiles far 
downstream, of two-stage turbine. 
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Figure 18.— Static pressure for two-stage 
turbine. 


on annulus wall boundary layers shows that it should be possible to 
improve the predictions in the region of the end walls. 

The turbine was fitted with static pressure tappings in the annulus 
walls. Figure 18 shows comparisons of observed pressure distributions 
with the theoretical calculations. The static pressure ratio is defined as the 
ratio of local static pressure to turbine inlet stagnation pressure. The main 
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point to note is the presence of an inverse pressure gradient in the region 
of the second-stage stator blade — static pressure on the inner wall greater 
than that at the outer wall — which was successfully reproduced by the 
through-flow analysis. An alternative to the through-flow method is what 
is known as the streamline curvature duct flow method (ref. 5). Frost 
(.ref. 24) has found that this method, which is widely used throughout the 
aircraft industry, did not predict the inverse pressure gradient. This 
example serves to show that when calculating the detailed internal 
aerodynamics, the flow inside the blade rows must be analyzed if a fairly 
accurate solution of the flow pattern is required. 


Single-Stage Axial Flow Turbine 

This single-stage, lightly loaded turbine was designed and tested at 
NGTE (ref. 25). In the initial through-flow analysis, no allowance was 
made for annulus wall boundary layers and the local polytropic efficiencies 
were assumed to be constant throughout the flow field. The predicted 
velocities (fig. 19) at turbine exit were in fair agreement with the observed 
values. Some measure of improvement in the region of the outer annulus 
wall was obtained by Herbert et al. (ref. 26) by allowing for the blockage 
caused by the boundary layers on the annulus walls. Improved matching 
of the experimental and predicted velocity profile would require a detailed 
boundary-layer analysis along the lines suggested by Gregory-Smith. 


EXPERIMENT 



MATRIX 

THROUGH 

FLOW 


Figure 19. — Axial velocity profile jar 
downstream of single-stage turbine. 



62 THEORETICAL PREDICTION OF FLOWS IN TURBOMACHINERT 

Low Pressure Ratio Centrifugal Compressor 

This centrifugal compressor was designed and tested by a firm in the 
United Kingdom. The results shown in figure 20 are the experimental and 
theoretical distributions of static pressure ratio along the shroud. The 
static pressure ratio is defined as the ratio of local to inlet static pressure. 
In performing the initial calculations, the values of local polytropic 
efficiency were assumed to be constant throughout the flow field and the 
slip factor equal to unity. The solution, although giving the correct trend, 
is in poor agreement with the observed pressures. By assuming a non- 
uniform distribution of local poly tropic efficiency and a slip factor of 0.91, 
the matching between experiment and theory was improved. This example 
shows that if a scientifically based model for the flow loss can be formu- 
lated then the through-flow theory might eventually be used to provide a 
quantitative picture of the flow pattern. 

High Pressure Ratio Centrifugal Compressor 

This example of a centrifugal compressor has been included to illustrate 
the use of the through-flow program at the design stage of a machine. The 
initial and modified (final) hub-shroud profiles are shown in figure 21. 
The only difference between the two impellers is that for the modified 
machine, the inducer extends beyond the leading edge of the splitter 
vanes, thus giving a deeper inducer section. The relative Mach number 
distributions along the hub and shroud profiles are shown in figure 22. It 
will be seen that the severe velocity gradient in the region of the inducer 


Figure 20 . — Sialic pressure distribution 
along the shroud of a low pressure 
ratio centrifugal compressor. 
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Figure 21 . — Hub-shroud profile of high pressure 
ratio centrifugal compressor. 
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Figure 22. — Hub-shroud Mach num- 
ber distributions for high pressure 
ratio centrifugal compressor. 


leading edge for the original impeller is to some extent alleviated in the 
modified impeller. This example demonstrates the effects of modifications 
that are possible within the limits of the same inlet and outlet areas and 
overall length of the machine. A new design may, of course, permit 
variations on all these factors and the use of a computer method helps in 
choosing the best combination. 

CONCLUSIONS 

Computer solutions for the meridional and blade-to-blade flow patterns 
in turbomachines have been described. The theory is based on the earlier 
work of Wu (ref. 11) and the numerical solution is obtained by finite 
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difference approximations to the governing equations. The main con- 
clusions are the following. 

Blade-to-Blade Flow 

(1) Comparisons with exact cascade solutions show that the blade- 
to-blade program gives an accurate solution for incompressible flow. 

(2) Analysis of turbine rotor blade sections shows that detailed flow 
patterns can be obtained which would probably be very difficult to find 
experimentally. 

(3) A comparison with experimental data for a turbine stator blade 
shows that the method gives a good estimate of high subsonic flow. This 
analysis demonstrates that the asymmetric finite difference grid developed 
here is an advancement over the conventional square or rectangular grid. 

(4) An example of a two-stage turbine illustrates that the three- 
dimensional pressure distributions can be predicted quite well. 

Meridional Flow 

(1) The matrix through-flow theory has enabled significant advances 
to be made in calculating meridional flow patterns. An analysis of a two- 
stage turbine shows that the theory gives a good estimate of annulus wall 
static pressure distributions. 

(2) An example of a centrifugal compressor shows that small modifica- 
tions to the impeller can have significant effects on the flow field. This 
analysis demonstrates that computer methods can help in selecting the 
“best” geometry. 

(3) A simple calculation of annulus wall boundary layers for a single- 
stage turbine enables the through-flow predictions to be improved by 
allowing for the blockage caused by the boundary layers. 

(4) Improved matching between experimental and predicted flow 
profiles depends on finding a better loss model and an accurate solution 
for the boundary-layer development along the annulus walls. 
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LIST OF SYMBOLS 


Q Mass flow 

Radial, axial, and circumferential coordinates 
V Absolute velocity 

TXT + 

** ivciauvc vciuuiuj 

x, y Coordinates with tilted axes 
X, a, y Flow angles 
^ Stream function 

Subscripts 

d Far downstream of blade row 
u Far upstream of blade row 

x ^-component 

y ^-component 

d> Circumferential component 
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DISCUSSION 


T. KATSANIS (NASA Lewis Research Center) : Mr. Smith has shown 
that approximate three-dimensional solutions for flow through a turbo- 
machine can be obtained by a combination of meridional and blade-to- 
blade solutions. This is fairly well known. However, we have here a fair 
number of examples showing both the strengths and weaknesses of these 
methods in applications. 

Limitations of the method should be clearly stated. It appears that the 
flow must be absolutely irrotational, steady relative to the rotating blades, 
and nonviscous, and that the flow must be complete subsonic. There must 
be other assumptions. Certain advantages are stated for the method and 
the program, but the basic assumptions and limitations are not specified. 

The author should specifically state how his method differs from those 
currently available. For example, a nonorthogonal mesh is used, but the 
corresponding finite difference equation is not given. Another example is 
equation (8), where a term has been added, but no explanation of its 
significance or why q(x,q>) is not also a function of l)\p/dx. 

Some comments must be made on one error. This is the statement that 
the boundary condition of zero velocity normal to the blade is not satisfied 
in references 8 and 9. This statement is not true. Further, the statement 
is made that finite difference stars with short limbs leads to large trunca- 
tion errors. It is true that the standard finite difference equation for 
unequal spacing has a larger truncation error than with equal spacing. 
However, this does not mean that the error in the solution will be larger. 
In fact, with a rectangular mesh, there is theoretically no loss in the 
accuracy of the solution due to an irregular boundary. This has been 
amply demonstrated by extensive use and experimentation with the 
programs of references 8 and 9. 

Mr. Smith does not discuss the reason for the use of iterative or relaxa- 
tion methods for solving matrix equations. There are two main reasons, 
one being the numerical stability which can be controlled by using a 
suitable rigorously calculated overrelaxation factor, which assures 
numerical stability with an optimum rate of convergence. The other 
reason is economy of storage, which is not shared by most direct methods. 
Certainly numerical stability cannot be improved with a direct method. 
The advantage of a direct method would be in reducing the computer 
calculation time. This reduction in computer time could conceivably be a 
real advance, provided that storage requirements are not significantly 
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increased, and that, as indicated, numerical stability is achieved. I would 
be very much interested in seeing a comparison of computer times for Mr. 
Smith’s program with the times for a method using optimum over- 
relaxation. 

M. D. WOOD (Cambridge University) : The paper indicates the 
magnitude of the recent advances made in calculating fluid flow in turbo- 
machines. The examples given show that compressible flow in turbine 
configurations can be predicted to a high degree of accuracy. However, 
as soon as compressor-type machinery is considered, the position is not so 
satisfactory. No one is really surprised, because the influence of boundary- 
layer growth and separation in compressor flow is likely to introduce 
effects which are of dominant importance. These effects are only recog- 
nized in the Wu equations through the presence of losses, and in general 
even these losses represent average or “smeared” values taken over 
appropriate computing planes. 

It is clear that there are few shortcomings in the equations of motion 
which Wu manipulates — the shortcomings are only in the simplifications 
we impose in order to obtain quick gains in current predictive accuracy. 

I therefore suggest that we should now have the courage to involve our- 
selves in combining the currently developing detailed calculations of the 
viscous effects in turbomachinery with the type of basic Wu program 
described by Mr. Smith. To take examples, we can see how boundary- 
layer separation in blade corners will lead to warping of stream surfaces. 
This warping can, in principle, be incorporated in the Smith-type pro- 
grams. Again, incorporating the predicted development of the boundary 
layer on blade surfaces would give better understanding of the “slip” 
factor for use in investigations of centrifugal compressors. Finally, in- 
clusion of the turbulent diffusion effects between adjacent fluid layers 
would lead to more realistic representation of the fluid forces in the Wu- 
type equations. 

Although this sounds like a daunting program, it is no more daunting 
than the thought, 10 years ago, of putting Wu’s equations on a computer. 
Perhaps the author would put my hopes into perspective by explaining 
what he intends to do next. 

A. S. MUJUMDAR (Carrier Corporation) : As pointed out by Dr. 
Katsanis, since the program itself is apparently the major contribution 
made by the author, it is unfortunate that it cannot be released for 
publication. Any comparison with the generally available Katsanis 
programs must, therefore, remain one-sided. The overrelaxation pro- 
cedure using proper grid spacing and an optimized overrelaxation factor 
should yield numerically stable results for well-guided geometries. As 
suggested by Wilkinson (ref. D-l) , the maximum relative velocity change 
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between successive iterations should be taken as the criterion for con- 
vergence rather than the maximum streamline deviation chosen by 
Katsanis. 

Since reference 12 in the author’s paper is not readily accessible, may I 
suggest that the finite difference analysis and the numerical scheme be 
included as an appendix to the paper when it is published. To my knowl- 
edge there is no “conventional” finite difference scheme to solve the 
Poisson-type partial differential equations; a number of variations are 
possible. 

Referring to figure 15 of the paper, could the author explain why the 
two-dimensional matrix solution appears to give better agreement than 
the three-dimensional solution with the experimental data for the pres- 
sure surface. 

Finally, I wish to bring to the author’s attention the experimental 
study of the flow in the blade passages of a radial turbine reported by 
Glenny (ref. D-2), which may be used to provide further checks for the 
computer code. 

R. C. DEAN (Creare Inc.) : I’m a little bit disturbed by perhaps the 
implications that you suggest about the use of potential analysis in 
centrifugal compressors. In my experience, the potential analysis usually 
considerably overpredicts the pressure rise in the wheel and underpredicts 
or predicts a low relative Mach number at the discharge of the impeller. 
We have found this through several comparisons between these solutions 
and data. The potential analysis you are suggesting implies, I think, that 
the flow follows the blading. I think it is very misleading to think that 
such a solution would work toward the back of the impeller. The important 
physics of the flow are not included in the analysis. 

L. MEYERHOFF (Eastern Research Group) : I have a number of 
questions. 

(1) Do you have any convergence criteria? 

(2) Are you able to predict the number of iterations for the con- 
vergence criteria? 

(3) How do you determine the trailing edge flow angle? 

(4) Did the addition of the E term referred to in your paper still keep 
the equations set up by Wu, exact? It is not clear whether the equations 
are still exact after you add these E terms. 

(5) Do you know of any analytical proof of the truncation error for 
the overrelaxation referred to in your paper? 

SMITH (Author) : The author thanks the five discussors for their 
review of this paper. 

First, taking the specific points raised by Dr. Katsanis, the purpose 
of this paper was to present numerical solutions for several turbomachines 
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to indicate the magnitude of the recent advances made in calculating the 
fluid mechanics rather than a detailed account of the mathematics of the 
flow models. The limitations of the methods and the finite difference 
approximations have been published in references 4 and 10. 

On the question of boundary conditions, I have perhaps not made my 
point clear. In Dr. Katsanis’ method (refs. 8 and 9) the flow domain is 
covered with a rectangular or square grid (fig. D-l) . To obtain the blade 
surface velocity at points such as A, the circumferential component of 
velocity, W is obtained from the relationship 


W+= — 


1 

bp dm 


and the resultant surface velocity W is determined so that there is zero 
velocity normal to the blade; thus, 



cos 0 


(D-l) 


where 0 is the local blade surface angle. For points such as B the meri- 
dional component of velocity, W m , is obtained from the relationship 


L^t 

bpr d<t> 


W m =—-r: 


and the resultant velocity is given by 


W= 


sin 0 


(D-2) 


I agree with Dr. Katsanis that these relationships ensure zero velocity 
normal to the blade. However, a restriction is placed on the resultant 
velocities; equation (D-l) is limited to | 0 |<60° and equation (D-2) is 
limited to | 0 [ >30°. This, I feel, implies an error in the derivatives of the 
stream function or the stream function values and when the condition 
of zero normal is imposed gives rise to an error in the resultant velocity. 



BLADE 

SURFACE 


Figure D-l. — Square finite difference grid. 
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Turning to the problem of the solution of the matrix equations, I must 
admit that I have no experience of iterative or relaxation methods. The 
method I have adopted is a direct method developed by the National 
Physical Laboratory (NPL), England. This approach, which removes one 
possible source of divergence, has proved to be very stable numerically 
and. as T pointed out, in my paper, only the band of nonzero elements are 
formed and stored in the computer. In the recent blade-to-blade com- 
puter program, the storage requirements have been further reduced by 
making use of magnetic tapes as backing store; if the width of the band 
matrix is W, the core store required is approximately 8JT 2 . There is a 
further point concerned with the basic techniques of the method of 
solution. Unlike “conventional” direct methods, the NPL procedures do 
not “invert” the matrix on every iteration; the solution is obtained by a 
backward and forward substitution process. The matrix equation to be 
solved is 

m-M-Cfl] 

The first step is a decomposition of the matrix [M~\\ thus 

[L]-[U] •[>] = [<?] 

where the matrices [ L ] and [t/J are lower and upper triangular band 
matrices, which are only computed on the first iteration. The solution is 
then obtained by (1) a process of forward substitution, solving for [Z] 
from 

[LHZ]=[ ? ] 

and (2) a process of backward substitution, solving for from 

[U]-M=[zj 

The direct method provides an exact solution for the matrix equations 
and so it could be argued that, since the overall process for finding the 
stream function distribution is an iterative procedure, it is not necessary 
to obtain an exact solution on the earlier iterations. It may well be that 
the best approach is a relaxation method on the earlier iterations, making 
no attempt to reduce the residuals to zero on each iteration, followed by a 
direct method on the final iterations. 

In answer to Dr. Wood, I would agree that we should now have the 
courage to extend the type of calculations I have described to include 
viscous effects. However, as is inevitable with an advanced calculation 
procedure, my experience of the use of the matrix methods has shown 
that, for the computer programs to become basic design tools, effort is also 
required in generating supporting programs for preparing geometric input 
data and graphical display of output data. This is one aspect I intend 
to examine. 
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I am grateful to Mr. Mujumdar for drawing my attention to experi- 
mental investigation of a radial turbine. With regard to the solutions for 
the turbine stator (fig. 15) I am unable to provide an explanation for the 
two-dimensional solution being in better agreement with the pressure 
surface experimental data than the three-dimensional solution. 

I agree with Mr. Dean that, in the case of centrifugal compressors, some 
of the important physics of the flow are not included in the analysis. 
However, I feel that, even with the assumption that the flow follows the 
blading, the computer tools can help the designer in selecting the best 
geometry. If a boundary-layer analysis of a potential flow velocity dis- 
tribution indicates, for example, separation in the inducer of a centrifugal 
compressor then I am sure Mr. Dean would agree with me that the 
designer would modify the geometry to overcome this problem. A number 
of examples illustrating this point are given by Dallenbach (ref. D-3) 
and Ball et al. (ref. D-4) . 

In answer to Mr. Meyerhoff, it is difficult to establish a unique con- 
vergence criterion. The criteria I have adopted are 

TOL = - — (through-flow method ) 


TOL = 


'J'—'l'i.-i 

Q 


(blade-to-blade method) 


It has been found that TOL can be reduced to 0.001 in 15 iterations for 
the through-flow method and 0.0001 in 10 to 30 iterations for the blade- 
to-blade method. 

The calculation of the flow angle is a problem I have avoided by as- 
suming it can be determined from existing empirical rules for the devia- 
tion. Clearly this is unsatisfactory, as I have indicated by the turbine 
example of figure 10. In a real flow the circulation is determined by viscous 
effects, particularly for compressor-type machinery as illustrated by 
figure 16. This is one aspect of turbomachinery fluid mechanics that 
demands research. 

The addition of the E(dip/dx) does not change the basic Wu equations. 
This term was added to keep the width of the band matrix to a minimum. 

Finally, on the question of the truncation error for conventional finite 
difference analysis, which Dr. Katsanis also raised, I can perhaps best 
illustrate my point by considering a square grid. For such a grid a five 
point star (refs. 8 and 9) is adopted to represent the Laplacian operator. 
Consider a star near to a boundary (fig. D-2) with one irregular limb. It 
will be supposed that for the star center, point 3, 



= 2 Unfn+T 
n= 1 


(D-3) 
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Figure D-2. — Five-point star with one 
irregular limb. 


where T is a truncation term and a u a 2 are the weightings. The function/ 
can be expanded as a two-dimensional Taylor series about the point 3 at 
which the function has the value f 3 and derivatives fz, x ,fz. v , fz.xx, etc. 
Substituting the Taylor series expressions into equation (D-3), it follows 
that 

/ h 2 h 3 \ 

fz,xx~\~fz,w ^1 (/a kfz,x~\~~ fz.xx T“/8,**x“t" .... I 


(fz~) 


h? h 3 

+<12 [fz~ hfz, v +—f3. n — — /s.rov-H ■ 


+03/3 


(f. 


k 2 k 8 

+ «4 (/3 + fe/3,v + “/3,l/l/ + “/3 .1/1/11 d - • 


/ h 2 h 3 

+ «5 (/3 + ^/3,x+“/3,xx + “/3.x*x"f‘ • 


-) 

-) 

) 


+r 

Since / is a general function, it follows that the coefficients of / and each 
of its derivatives may be equated on each side of the above equation. 
There are five disposable constants a,(i= 1(1)5) and so only the coeffi- 
cients of /, f x , f v , f xx and f V y may be used in order to make the finite differ- 
ence approximation independent of the low-order derivatives. From the 
coefficient equations, it is easy to show that the solution for the weightings 
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2 

(7 + l)/t 2 


a 3 = 


2(y + 1) 

7 h? 


&4 " 


(7+l)7/i 2 


a s - 


1 

"h* 


where 7 = k/h. 

The truncation error is 

h(y 2 —l) 

T'= f3,yyy-\-0(h%,yyyy) (D-4) 

3(7 + 1) 

It is seen, therefore, that as the irregular limb gets shorter (i.e., 7 de- 
creases) the truncation error increases. 
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